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ABSTRACT 

We introduce a differential equation for star formation in galaxies that incorporates 
negative feedback with a delay. When the feedback is instantaneous, solutions approach 
a self-limiting equilibrium state. When there is a delay, even though the feedback is 
negative, the solutions can exhibit cyclic and episodic solutions. We find that peri- 
odic or episodic star formation only occurs when two conditions are satisfied. Firstly 
the delay timescale must exceed a cloud consumption tiniescale. Secondly the feed- 
back must be strong. This statement is quantitatively equivalent to requiring that 
the timescale to approach equilibrium be greater than approximately twice the cloud 
consumption timescale. The period of oscillations predicted is approximately 4 times 
the delay timescale. The amplitude of the oscillations increases with both feedback 
strength and delay time. 

We discuss applications of the delay differential equation (DDE) model to star 
formation in galaxies using the cloud density as a variable. The DDE model is most 
applicable to systems that recycle gas and only slowly remove gas from the system. 
We propose likely delay mechanisms based on the requirement that the delay time is 
related to the observationally estimated time between episodic events. The proposed 
delay timescale accounting for episodic star formation in galaxy centers on periods 
similar to P '--^ 10 Myrs, irregular galaxies with P ^ 100 Myrs, and the Milky Way 
disk with P ^ 2 Gyr, could be that for exciting turbulence following creation of 
massive stars, that for gas pushed into the halo to return and interact with the disk 
and that for spiral density wave evolution, respectively. 



1 INTRODUCTION 



Gas present in a galaxy fuels star formation or nuclear 
black hole growth. However both star formation and ac- 
tive galactic nuclei then release energy and momentum 
into the interstellar medium (ISM). Consequently the ac- 
tivity can suppress subsequent star formation. The pro- 
cess in which part of the output of a system is returned 
to its input and influences its further output is termed 
"feedback." Early studies showed that when feedback by 
radiative heating is taken into account during gas accre- 
tion onto a central mass, steady solutions may not exist 
l|Ostriker et al.l 1 19761 ) and the feedback proce ss can cause 
oscill ations or periodic bursts of accretion (|Cowie et al.l 
1 197^ ). Simulations taking into account feedback processes 
illustrate that gas flows and star for mation in galaxies 
can exhibit episodic or cyclic behavior (iDong et al.l [200^: 
Pelupessv et al.ll2004 ICiotti fc Ostrikerll2007l: IStinson etal] 
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oTphin et al.l l2003l : 
Dellenbusc h et al.l 

disk ( Rocha-Pi nto et al.l 

llBland-Hawthorn fc CohenI l2003l : IWalcher et all l2006l : 
ICecil et all I2OOII) and the statistics of distant galaxies 
i Glazebrook et al] Il999l ) infer that multiple events of 
vigorous star formation, separated by millions to billions 
of years, can occ ur even in i s olated galactic systems. 
Other studies (e.g.. Ivan "Ze3l200ll : [s"killmanll2005^ find little 
evidence for episodic star formation. However, theoretical 
work has primarily fo c used o n self - regul a ted star formation 
jAnderscn & Burkcrt' '2000'; 'sil^ 2001; 'Elmcgrccn' '200^; 
Monaco 2004; Krumholz ct al. 2006; Slyz ct al. 
Li et al.l l2006l: iDib et all [2O O6: Jou ng fc Mac Low! 
Elmegreenll2007(: iBooth et al.li2007.: .Wada fc NormanI 



Schave fc Dalla Vecchial 20071 : Robertson fc Kravtsovl 



iUP' 

20071 ) or alternatively can asy mptotically converge onto a 



self-regulated equilibrium st ate (lAndersen fc Burkertll200ol : 
iRobertson fc Kravtsovll2008l ). 

Galaxies display complex star formation histories. 
Studies of irregular galaxy populations (e.g., iTosi et akl 



and has not explored when episodic rather than a steady 
rate of star formation is expected. 

As gas flows involving energy input, heating and cool- 
ing are complex, there is no simple way to predict when 
behavior is episodic or cyclic. However it is possible that 
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average quantities can be estimated for these flows and 
relations based on these quantities can be used to clas- 
sify their behavior. Delay differential equations can ex- 
hibit solutions that asymptotically approach a self-limiting 
equilibrium state and those that are periodic, even when 
feedback is negative. Consequently these equations can 
be used to differentiate between these two behaviors. De- 
lay differential equations have been used to model bi- 
ological systems with delayed negative feedback (e.g. , 
Wazewska-Czvzewska fc Lasota Il988l: Gvori fc Ladaslll99ll : 



Gurnev et al.lll980l : lKulenovic et al.ll 19891 ) but have not been 
applied to astrophysical systems. In this paper, using a delay 
differential equation, we determine when cyclic or periodic 
behavior is exhibited by the solutions rather than a smooth 
decay to a self-regulated steady state. We apply the theory 
to star forming galaxies, identifying delay mechanisms that 
could account for episodic accretion events inferred from ob- 
servations. 



2 ONE DIMENSIONAL FEEDBACK MODELS 

We begin by considering a galactic disk model with cloud 
surface density E(f) in units of mass per unit area that 
depends on time, t. This density could represent the total 
disk gas, or the gas in self gravitating clouds, or the gas in 
molecular form, depending upon the setting. The gas den- 
sity available for star formation decreases when clouds are 
dispersed following star formation. Conversely, the gas den- 
sity increases during accretion, coagulation or cooling, all of 
which can enhance star formation. We therefore write 

E(t) =p(S,t)-/i(E,i) 

where E — dT,/dt. Here the function gi(E,t) is the accre- 
tion or cl oud for mation rate. A Schmidt star formation law 
l|Schmidtt 1 959; Kennicuttl 19981 ') relates the cloud or gas con- 
sumption rate to the disk density variable with the function 
/i(E). In principle, the accretion rate also depends on time 
in a non-trivial manner. For example, it could depend on 
the previous star formation rate. We do not expect feed- 
back to be instantaneous as it takes millions of years for a 
burst of star formation to produce type II supernovae, and 
winds and supernova remnants require time to evacuate gas 
or induce turbulence in a gas disk. Following a burst of star 
formation, accretion onto the disk would not resume until 
heated, evacuated or dispersed gas has had time to cool and 
reform into clouds. 

Before introducing complicated functions for the accre- 
tion rate, we first consider the simplest case, that lacking any 
feedback, g{x,t) = A, corresponding to a constant accretion 
or cloud formation rate. The above differential equation can 
be written 



i = f{x) ^ A~ Bx" 



(1) 



where we have replaced E with the vari able x, use a Schmidt 
type star formation law (|Schmidtll 19591 ) with positive power 
index a, and A and B are positive constants. By setting 
dx/dt = and solving for x we find a fixed point, x,, 
corresponding to the self-regulated steady state or equilib- 
rium value at Xt — {A/B)^^". We can assess the nature 
of solutions by taking the derivative of the right hand side 



with respect to x; oi ^ = —Bax°'~^ . This derivative is al- 
ways negative and is smoothly decreasing function implying 
that solutions always smoothly (asymptotically) approach 
the equilibrium state solution on a timescale determined by 
the inverse of this derivative. There are no oscillating or di- 
vergent solutions. 



2.1 Instantaneous Feedback 

The case of instantaneous feedback can be modeled with the 
assumption that the accretion rate is affected by the current 
star formation rate, which in turn is set by the density of 
the disk. We expect that feedback would occur by reducing 
the quantity of gas available for star formation in the disk 
when the star formation rate is high. Since the gas quan- 
tity available to form stars is reduced by the star formation 
process, the feedback is negative. We can describe this sit- 
uation with an accretion rate g[x) = AG{x), where G{x) is 
function that approaches unity when x is small and there is 
no feedback, and drops to zero when x is large, star forma- 
tion is vigorous and the energy arising from it has prevented 
further accretion or cloud formation. A simple form for the 
function G that satisfies our requirements is G{x) — e~^^'~' 
for which C > 0. The parameter C depends on the star for- 
mation rate that is effective at cutting off accretion or cloud 
formation. 

The evolution of the disk density is then described by 

e/C 



X — f{x) — Ae 



Bx 



The equilibrium state can be found by solving x 
and satisfies 



a ^ - 

x^=^e 



(2) 
for X 



(3) 



The derivative of / 
df 



dx 



-Bax 



Since A,B,C,a > the derivative is always negative and 
solutions always smoothly approach the equilibrium state. 
Because the feedback is negative there is no instability, and 
no periodic or cyclic solutions exist. Solutions to this equa- 
tion resemble those that do not oscillate shown in in Figure 

m 

It is useful to define two timescales, a consumption 
timescale dependent only on the second term of equation 
[2] and evaluated at Xt 

" ' (4) 



dx 



BxT 



and the timescale to approach equilibrium, teq, that depends 
on the derivative of / evaluated at a;*. 
-1 



= icon {a + x,/Cy 



(5) 



It is also useful to quantify the strength of the feedback 
near the equilibrium point by looking at the sensitivity of 
the accretion term g{x) or 



dg X 



dx g(x) 



Xt 



(6) 



The above "strength parameter" is large when the feedback 
is strong. 



2.2 Delayed feedback 

When feedback is delayed the current accretion rate is re- 
duced by the star formation rate at an earlier time t — t 
where r is the delay timescale. The accretion rate is g{x{t — 
r)) and the model described by equation [2] becomes 

= fix, t) = ^e"'*-")/'^ - Bx{ty. (7) 

Terms of this form were considered by iDong et al. I l|2003l ). 
In the limit of r — > we recover equation [2] for in- 
stantaneous feedback. The above differential equation be- 
longs to the class of one dimensional equations with de- 
layed negative feedback which in cludes the delayed logis- 
tic equation, the model used by i Gurne y et al. (1980) to 
describe the dynamics of Nicholson's blowflies and the 
Lasota-Wazewska model for the survi val of red blood cells 
ijWazewska-Czvzewska fc Lasotalll988h . 

Even though the feedback is negative, the above differ- 
ential equation has oscillating solutions but not for all values 
of the 4 positive parameters A,B,C,t and index a. Figure 
[T] shows two solutions that converge to a periodic solution 
that oscillates about the equilibrium state forever and one 
that oscillates while decaying to the equilibrium state. The 
equilibrium state for this differential equation is identical to 
that for the equivalent model lacking delay and is a solution 
of equation |31 To display solutions we directly integrated 
Equation [7] with a first order or Euler method. We allow 
the delayed feedback to initiate only at times t > ta + t for 
initial time Iq. 

For non-extreme values of initial conditions and param- 
eters, the solutions exhibit 3 types of behavior: 

(i) The solutions lack oscillations. After some time period, 
solutions smoothly or asymptotically approach the stable 
equilibrium state. 

(ii) The solutions exhibit oscillations about an equilib- 
rium state but asymptotically approach that state. 

(iii) The solutions oscillate and are attracted to a periodic 
function or cycle. 

When oscillating solutions are present, the oscillation 
period is approximately four times the delay timescale or 
P ~ 4r. Roughly speaking, this follows by considering the 
equation x{t) = x{t + n/2) that has the solution x{t) = 
sin(a::) with a period of 27r. From figure[T]we see that the ac- 
tual period displayed by the oscillating solutions is approxi- 
mately 5r. This is broadly consistent with the approximate 
estimate for the period of 4r. 

To facilitate classification of solutions, we transform 
equation [7] into dimensionless form. We let a dimensionless 
density variable y = x/xt and time variable T = t/tcon with 
the depletion or consumption timescale defined in equation 
m Using these new variables, equation [7] becomes 

^ = e{i-!/(T-^))s_ . (8) 

where the dimensionless parameters are 

f = S = x./C, (9) 

and we have used equation |6] for the feedback strength S. 
For a given exponent, a, equation [S] only depends on two 
parameters, r and S, thus solutions of this equation can be 
classified based on estimates for these two ratios alone. 
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Figure 1. We show the results of integrating equation [7] for 
different dimensionless ratios S = x./C and f = r/tcon- The 
power index is o = 1.4. The solution with high values of these 
f and S is strongly non-sinusoidal and approaches a high ampli- 
tude periodic solution (thick solid line). At lower values of these 
parameters a lower amplitude but periodic solution is approached 
(thick dashed line). For even lower values there is an oscillating 
solution that decays to the equilibrium value (solid thin line). 
At the lowest values of f and S the solution asymptotically ap- 
proaches the equilibrium value (dotted thin line and dot dashed 
line). The ratios are listed in the plot key. The oscillation period 
is approximately five times the delay timescale. The j/-axis is x di- 
vided by that of equilibrium state, . The i-axis is time divided 
by the delay time, r. 

In the case of power index a = 1, the above DDE 
(equation [7| is the same as the model used to describe sur- 
vival of red blood cells by IWazewska-Czvzewska fc Lasotal 
(1988). Initially positive solutions of the dimensionless ver- 
sion (equation [S]) oscillate about the equilibrium value, 
= 1, if and only if 

r5'e("+'> > 1, (10) 

jKulenovic fc Ladaj[l987l : iGvori fc Ladaslll99ll ). The equi- 
librium value is a global attractor (solutions approach this 
value) when 

5(1-6"") < In 2 (11) 

jKulenovic et al.lll989l : iGvori fc La"da^ll991^ . If this condi- 
tion is not satisfied, a periodic oscillating attractor may ex- 
ist. 

A more generalized oscillation criterion that can be used 
when Q 7^ 1 is 

Ml = f Se<""+') > 1, 

where we have defined Mi as a parameter describing the 
nature of the DDE. We have derived this criterion in 
the appendix u s ing th e procedures rigorously described by 
ICvori fc Ladaj (|l99ll ) and following the example given for 
the Lasota-Wazewska model. By analogy with equation 1111 
we guess that j/* is a global attractor when 

Afa = 5(1 - e~"") < ln2. 

The parameters Mi , M2 > 1 when 

r = r/tcon >1 and 5 > 1. (12) 
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Figure 2. Amplitude contours are shown for solutions of equa- 
tion [7] after the solution has decayed either to an equilibrium 
value or a periodic cycle. The amplitude is the maximum divided 
by the minimum value of x during the cycle. The lowest con- 
tour divides solutions that converge to a periodic function from 
those that decay to an equilibrium value. This contour has the 
amplitude value 1.05. The remaining contours are evenly spaced 
in the log with amplitude values 3.16, 10.0, 31.6, and 100 (second 
from bottom to top contour). For the solutions to exhibit peri- 
odic behavior we find that the ratio of the delay to consumption 
timescales r/tcon ^ 1 and feedback strength S = x^/C > 1. The 
amplitude of the cycles increases with increasing S and r. a) For 
index a = 1.5 b) For index a = 1.0 



Using equation [5] we find that the second of these condi- 
tions is equivalent to a constraint on the timescale to reach 
equilibrium 

Thus when two conditions are satisfied we predict periodic 
solutions: 

(i) The delay timescale exceeds the consumption 
timescale. 

(ii) Feedback is strong and eflFective at shutting ofT accre- 
tion or cloud formation near the equilibrium density. Equiv- 
alently the timescale to approach equilibrium exceeds the 
consumption timescale by a factor similar to 2. 



We consider how the amplitude of oscillations depends 
on the parameters. We describe the amplitude of oscillations 
as the ratio of the maximum x divided by the minimum in an 
oscillation period, after the system has converged to a cycle. 
Oscillation amplitudes are shown in Figure [2] as a function 
of r and strength S, for index oi = 1 and a = 1.5. The 
further away from the line dividing asymptotically decaying 
solutions from those with periodic solutions, the larger the 
oscillations about the equilibrium value. The amplitude of 
the oscillations does not depend on the initial conditions but 
rather on the parameters defining the differential equation. 
Since x cannot cross zero when large oscillations are present, 
the periodic solutions are less symmetric or less like sinusoids 
but exhibit spikes followed by longer periods of low periods 
of accretion when the amplitudes are high (see figure[T}. This 
follows as the accretion rate depends on the exponential of 
X so when x is high, it can take a long time for the system 
to recover from a previous episode of star formation. 



3 APPLICATIONS OF THE ONE 

DIMENSIONAL MODEL TO GALAXIES 

Our DDE model for delayed feedback is appropriate if the 
mean gas density averaged over long periods of time is 
nearly constant. This follows because the form we have 
for the accretion or cloud formation rate does not change, 
though it does depend on the past disk density. The DDE 
model is best applied to systems that recycle gas and 
only slowly remove gas from the system. Star formation 
laws illustrate t hat st ar formation is inefficient. For exam- 
ple, [kennicuttj (|l998l ) found that star formation rates in 
nearby galaxies could be described with E ~ eEJl, where 
Q. is the a ngular rotation rate and the efficiency is low, 
e ~ 0.017 (|Kennicuttlll998l '). This suggests that we should 
not adopt as our defining variable the total density in molec- 
ular and atomic gas but rather that in molecular clouds 
or self-gravitating clouds as adopted in explanations for 
the Schmidt-Kennicutt star formation law and ob servational 
studies of mole c ular gas in galaxies (iGao fc Solo mon 200j; 
'Wu ct al.' 200 51: iKrumho lz ct al' 2006': iBlitz fc Rosolowsk^ 
[2006, : .Robertson fc Kravtsov. ,2008'). In this case the DDE 
tracks cloud formation and cloud disruption following star 
formation. 

Molecular clouds are estimated to last ~ 2 — 
3 X lO^yrs and are disrupted following star formation 
l|Blitz et al.ll2007h . Theoretical work suggests that clouds 
disperse a fter a few times their fr ee fall or dynamical 
timescale l|Krumholz fc McKed [20051 ) so lifetimes of star 
forming clouds could b e shorter in denser environments 
l|Wada fc Normanll2007l ). such as circumnuclear disks. The 
depletion term in equation [7] has B = and index ex — 1, 
so the consumption timescale in our model is the mean cloud 
lifetime, icon ~ tc- 



3.1 Possible Delay Timescales 

For star formation to be maintained in a disk, clouds 
must constantly reform. Enhanced turbulence in the disk 
should reduce the star fo r mation rate (e.g., ISilkl 1200 ll : 
Ide Avillez. fc Breitschwerdtl 12004 iLi et al.l 120061 ) . Turbu- 
lence increases the disk thickness reducing the mean 
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density an d inc r easing the mean f re e fal l timescale 

(cf.,^ fsiii3 I2OOII: iKrumholz fc McKed l2005l: iDib et all 



Mac Low~^200' 



McKee fc Ostrikej l2007[ 



2006: Joung 

Robertson fc Kr avtsov^2008i ). The primary energy source for 



the turbulence is expected to be from supernovae, though 
differential rotation, gravitational and magneti c instabili- 
ties and stellar outflows could also pl a y a role (|Silk 200 ll: 



Kim et al.l l2003l: iQuillen et al 



2007l : lde Avillez fc Breitschwerdtll2007h . Thus a delay time 



20051: iPiontek fc Ostrikeij 



for feedback is the sum of the time for massive stars to move 
off the main sequence and produce supernovae (a few times 
lO^yr), the timescale for the supernova remnants to reach 
their maximum size (of order 10^ yr but depending on the 
ambient pressure and density), and the timesca le for them 
to be mixed into the disk (e.g., IPib et al]l2006t ). This last 
timescale is a turbulence mixing timescale, tmix ~ h/a, that 
depends on the gas disk thickness, h, and gas velocity dis- 
persion, a. The mixing timescale is similar to a few times 
10^ yrs in the solar neighborhood. Thus the delay time for 
a reduction in the rate of molecular cloud formation by tur- 
bulence in disks such as the Milky Way is a few times 10 
Myrs and dominated by the timescale for mixing and super- 
nova remnant expansion. (The timescale could be shorter if 
star formation is triggered by the rapid collapse of the evac- 
uated region (~ 2 Myr) shortly after the hot gas escapes 
the disk.) Both turbulent mixing timescales and supernova 
remnant expansion timescales should be longer in the out- 
skirts of galaxies and in irregular or dwarf galaxies where 
the densities and pressures are lower. In contrast, on the 
scales of circumstellar disks (~ 10 pc), mixing and super- 
nova remnant expansion timescales should be shorter due 
to the higher densities and pressures and larger velocity dis- 
persions. 

In spiral galaxies, molecular cloud formation occurs pri- 
marily in spiral arms so their formation is triggered on a 
timescale related to the spiral density wave pattern rather 
than on a timescale related to turbule nt mixing of super- 
nova remnants (e.g., lElmegreenI 120071) . A possible longer 
delay timescale is that for s piral density waves to evolve 
(e.g., IClarke fc Gittind I2OO6D . When the Toomre Q pa- 
rameter is greater than 1.5, spiral structure is suppressed. 
Here Q = where k is the epicyclic frequency and 

G the gravitational constant. Th e Q parameter is related 
to the gas fre efall timescale (e.g., 'McK ee fc Ostrikeill2007l : 
IKrumholz fc McKee 2005 ) and so its value can be discussed 
in terms of a self-regulated star formation model. Spiral den- 
sity waves are expected to grow o n a times cale of a few rota- 
tion periods (^S cllwoo d fc Carlbe re 1984; Vorobvov fc The"i3 
I2OO6I : IClarke fc GittinsI I2OO6I ). Star formation not only in- 
fluences the gaseous velocity dispersion but lowers the mean 
stellar velocity dispersion and increases the stellar mass den- 
sity. Hence the current strength of spiral structure (set by 
Q) may depend on the star formation rate a few galactic ro- 
tation periods ago. In this setting the cloud formation rate 
would be forced by spiral arms sweeping through the disk 
with an oscillation period dependent on the spiral pattern 
speed and amplitude dependent on the strength of spiral 
structure. This amplitude would be the quantity that expe- 
riences the delayed feedback. 

A third candidate for a delay timescale is that for ma- 
terial driven out of the disk to return and stir the disk. This 
could be influenced by a cooling timescale for hot and low 



density gas in the galactic halo. This timescale would be 
longer than the local disk turbulent mixing timescale and 
would be of order 10* — lO" yrs. It may be related to the 
100-200 Myr relaxation timescale exhibited by simulations 
11 de Avillez. fc Breitschwerdtll2004l ; Ijoung fc Mac Lowll2006l : 
IStinson et al.ll2007l ) but could also depend on the dark mat- 
ter halo mass or density (as discussed in these works). 

In summary, the relevant consumption timescale is the 
molecular cloud lifetime of order 10 Myrs but could be 
shorter in denser environments. For delay timescales we have 
three primary candidates: 1) The timescale for supernovae 
to enhance disk turbulence (a few times 10 Myrs but longer 
at lower densities and pressures). 2) The timescale for gas 
heated up and moved into the halo to cool back into and 
stir the disk (order 10® — 10^ yrs). 3) The timescale for spi- 
ral arms to evolve (a few times the rotation period). Fu- 
ture work may identify delay times associated with other 
processes such as magneto-gravitational instabilities, or in- 
ternally generated stellar outflows. The delay timescale as- 
sociated with disk turbulence may not exceed the cloud 
consumption timescale. However delay timescales associated 
with larger scale turbulence and cooling in the halo and spi- 
ral arm evolution are likely to exceed the cloud consumption 
timescale. 



3.2 Delay mechanisms as suggested by 
observations 

We now put these timescales in context with observations 
keeping in mind that the DDE (equation[7| displays episodic 
bursts only when the delay timescale is longer than the con- 
sumption timescale. 

The survey by iRocha-Pinto et al] (|2000al ) reveals that 
star formation in the solar neighborhood experienced 3 
bursts each separated by about 3 Gyrs. A delay timescale 
of one quarter of this or about 0.8 Gyr would be required 
to predict this periodicity with the DDE of equation [T] 
As spiral structure is responsible for molecular cloud for- 
mation in the solar neighborhood a possible delay mecha- 
nism is the timescale for spiral arms to evolve. The time 
0.8 Gyr corresponds to 3 rotation periods at the solar circle. 
IClarke fc Gittind (|2006l ) have previously proposed that vari- 
ations in spiral arm strength could affect the star formation 
rate. Here we couple the gas and stars, relying on feedback 
and a delay time but involving the same principle, that the 
spiral density waves are a strong trigger for star formation. 

Surveys of galaxy centers have revealed that most late 
type and elliptical galaxies harbor circumnuclear star clus- 
ters dBoker et al.l I2OO2I: iKoda et al. l l2005l : ICote et al.|[2006l : 
IChristopher et alj I2OO5I ) and have experienc ed star forma- 
tion in their nuclei in the past few to lOOMyr ([Veilleux et al.l 
'1994': 'Bland-Hawthorn fc Cohen"20p^ ; IWalcher et al ] |2006l : 
Quillen ct al. 2006; Cecil et al. 200l|). The sizes of these star 
clusters ranges from tens to a few hundred pc and gas den- 
sities of 10^-IO^Mq pc~'^. Since the gas densities are high, 
cloud lifetimes should be shorter than that for molecular 
clouds in the Milky Way's disk or Local Group galaxies. Su- 
pernova remnant expansion and turbulent mixing timescales 
may be shorter than in the solar neighborhood due to higher 
pressures. However the timescale for stars to evolve must be 
similar in both settings. We expect episodic star formation 
with a period similar to a few times 10^ years (set by stellar 
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evolution of massive stars). This behavior would only occur 
when the timescale for excitation of turbulence in the disk, 
depending on the timescale for stars to produce winds, is 
longer than the lifetime of the star forming self-gravitating 
clouds. 

Studies of irregular dwarf galaxies have revealed 
that they have complex star formation histories expe- 



riencing separated bursts of star formation separated 


by a hundred Myrs to Gyrs (e.g., 


Tosi et al. Il99ll: 


Dohm-Palmer et al.1 2002 


: DolDhin et al. 


I2OO3I: ISkiUmanI 


20051: YouuE et alj [20071: 


Dellenbusch et al.' '2008"). Recent 


simulations (Peluoessv et al. 2004: Stinson ct al. 200711 have 



illustrated periodic bursts of star formation separated by 
200-400 Myr. The simulations do not display strong spiral 
structure. The sp i ral str ucture mediated model proposed by 
IClarke fc GittinsI (|2006l ) can account for bursts of star for- 
mation in dwarf galaxies, however this model cannot ac- 
count for the bursts seen in these simulations as they lack 
spiral structure. The delay timescale must be one quarter 
of the time between bursts or 50-lOOMyr. The supernova 
remnant expa nsion timescale for the galaxy simulated by 
IPelupessv et a l. (2004) is similar to that of a supernova in 
the solar neighborhood as the interstellar medium pressures 
are similar. Likewise turbulent mixing timescales are sim- 
ilar. Hence the long inferred delay timescale must involve 
longer timescales such as for cooling of material in the ha- 
los of these galaxies and interactions between this cooling 
material and the disk. 

In all three of these cases, it is likely that the delay 
timescale exceeds the cloud consumption timescale, one of 
the conditions for the DDE to exhibit cyclic solutions. We 
base our choices for the likely delay mechanism on the re- 
quirement that the delay time is related to the observational 
inferred timescale between episodic events. Thus we suspect 
that the relevant delay timescale accounting for episodic star 
formation in galaxy centers, irregular galaxies and the Milky 
Way disk could be that for exciting turbulence following cre- 
ating of massive stars, that for gas pushed into the halo to 
return and interact with the disk and that for spiral den- 
sity wave evolution, respectively. In all three cases, the total 
supply of gas is consumed only slowly leaving a reservoir for 
ongoing star formation. Since the feedback is delayed on a 
timescale that exceeds the cloud consumption timescale, re- 
current and periodic star formation events could occur even 
though the feedback is negative. 



3.3 Is the feedback strong enough? 

We now discuss the second requirement for cyclic solutions, 
that feedback be effective at reducing the formation rate 
of molecular clouds. We have characterized the feedback 
strength, S, with a parameter defined in equation [6] that 
describes the change in cloud formation rate caused by a 
change in cloud density. Only when 5 > 1 are the solutions 
to the DDE periodic in behavior. Consequently we need to 
estimate the change in the cloud formation rate (or star 
formation rate) caused by a small change in the mean gas 
density. 

There are few references that have considered th e 
timescale for cloud formation (q.v. iPadoan et al.l |2006| ). 
More commonly, a density spectrum resulting from tur- 
bulence has been used to predict the number of clouds 



above a critical density. The star formation rate is esti- 
mated from this gas fra ction divided by the dynamica l 
timescale at that de nsity (Elmcgrcen 2002; 'Kra vtsov|[2003l : 
iKrumholz fc McKeQ,2005. : Wada fc Norman 200'^. A nearly 
universal property of isothermal turbulent media in exper- 
imental and numerical simulation studies is that t he cloud 
densities have a l og normal density distrib ution (jWarhaftl 
l2000l : |Pum ir"l994':'Padoan fc Nordlund'2002'). We adopt this 
distributiorQ to estimate the strength parameter S in equa- 
tion [6] 

Stars are born primarily in the densest clumps that form 
as a result of turbulence within the interstellar medium. 
The disk velocity dispersion is predicted t o be propo rtional 
to the square root of the supernova rate (|Dib et al.l 12006. ). 
So the mean gas density should depend on the square root 
of the star formation rate. The star formation rate is esti- 
mated from the fraction of mat erial in the densest clumps 
or that above a critical density. jP adoan fc Nordlund 2002 
Elmegreenll2002l: IKrumholz fc McKe6,2005. : ,Kravtsov,2003 



Wadafc Normanll2007f ). The fraction of the mass with a 



density, p, larger than a threshold, pc 



fc 



/p7 PPip)dp 



/o°° PPip)dp 

where the normalized probability density function 

dlnu 



p{u) = (27rA^)"^''%xp(-0.5[lnu-lnuo]VA^)' 



du 



(13) 



(14) 



Here u = p/p is the density in units of the mean density 
and ln«o is the mean of the normal distribution. The mean 
and dispersion of the normal distribution depend on the 
Mach number on the large st scale and are in the range 1-5 
(|Padoan fc Nordlund|[200a ). 

After integrating, we estimate oc erfc ( ^'"^^^--^ ) 

(based on equation 20 bv IKrumholz fc McKeeir2005h . where 
the critical density ratio Ucrit = Pc/p, we have used a com- 
plementary error function and assumed that the critical den- 
sity ratio exceeds the mean by more than a few disper- 
sion lengths A. In the large asymptotic limit this becomes 
fc ~ Q-O-'^^orit/'^) _ ^ change in the density ratio Ucrit leads 
to a change in the cloud fraction 



S- 



dfc u 
du fc 



2 In Ucrit 



(15) 



The above ratio, equivalent to the strength parameter de- 
fined in equation[6l tells us how large a change in the fraction 
of clouds above the critical density is caused by a fractional 

crit IS eS' 

timated to be in the range of 10* - 10'' 
Krumholz fc M cKcc 2005). For A = 2.4 



Elmegreen 


2002 


Elmegreen 


2002 



Padoan fc Nordlundii2002l ) and Ucrit = IQ" 



the above frac- 
tion 5 ~ 4. We expect the condition strength 5* > 1 for 
our model is satisfied but that the strength is also not ex- 
tremely large. For delay times exceeding the gas consump- 
tion timescale by a moderate factor with 5* ~ 4 we would 

^ While there is no theoretical basis for this distribution, R. 
Sutherland (personal communication, 2008) points out that it is 
a natural consequence of a turbulent cascade with multiplicative 
rather than additive random phases due to folding and stretching 
within the medium. 
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predict solutions with moderate amplitude oscillations (see 
Figure [2Jd). 

The feedback strength estimate shown in equation [15] 
suggests that the feedback would be weaker at higher Mach 
number but stronger at lower mean density, if the critical 
densit y is similar in different environments. IStinson et al.l 
(|2007l ) found that oscillations were lower amplitude for 
larger simulated dwarf galaxies. Figure [5] showing the ampli- 
tude as a function of feedback strength and delay timescale 
implies that the feedback strength would be lower for the 
larger simulated dwarfs because they have longer delay times 
and because their mean gas density is higher. 

Further examination of these simulations may test the 
hypothesis that equation [15] describes the feedback strength 
and is consistent with the relationship between oscillation 
amplitude and feedback strength predicted by the model. 
The above estimate for the feedback strength is indirect as 
we have used a steady-state star formation rate to estimate 
the cloud formation rate. Timescales displayed by simula- 
tions of th e density evolution and m olecular cloud forma- 
tion (e.g., iGlover fc Mac Low! 120071 ) might allow a better 
and more appropriate estimate for the feedback strength. 
The strength we estimate above was based on a local prob- 
ability density distribution but when feedback delay is very 
long (such as suggested in the solar neighborhood) the cloud 
formation rate should be integrated azimuthally around the 
galaxy and across spiral arms. 



4 SUMMARY AND CONCLUSION 

It is now widely recognized that a detailed understanding 
of feedback and accretion processes is essential to progress 
in many fields of astrophysics and across the entire cosmo- 
logical hierarchy, from galaxy clusters down to the scales of 
individual star forming regions. In order to progress, we will 
need huge improvements in analytic algorithms and com- 
puter power, as well as better conceptual tools for classifying 
complex behavior. Some processes may indeed be episodic 
or cyclic, while other instances may exhibit quasi-periodic 
cycles on the way to fully chaotic behavior. A deeper un- 
derstanding requires that we should to some degree be able 
to distinguish between these two very different dynamical 
manifestations for open and closed systems. 

Here we have introduced a simple differential equation 
model that captures some of the complexity exhibited by 
astrophysical star forming systems with feedback. We intro- 
duce a one dimensional DDE for the molecular cloud density 
that allows cloud formation to depend on the star formation 
rate but at a previous time. Thus current star formation only 
affects the cloud distribution at a future time, we denote the 
delay time. The feedback is negative, so in the absence of 
delay there are no cyclic solutions or instabilities and all 
solutions asymptotically approach a self-limiting value. 

We illustrate that even when the feedback is negative a 
delay can cause cyclic or episodic behavior. The DDE cap- 
tures phenomena exhibited by astrophysical simulations of 
this process, including periodic solutions in some cases but 
not in others. The DDE allows us for the first time to clas- 
sify the solutions and predict when an astrophysical system 
is self-limiting or likely to exhibit periodic behavior based 



on timescales that are related to physical feedback and star 
formation processes. 

We find that periodic behavior is likely when two con- 
ditions are met. First, the delay timescale must exceed the 
cloud consumption timescale. Secondly, the star formation 
must be effective at reducing the rate of formation at den- 
sities near the self-limiting or steady state value. This is 
equivalent to requiring strong feedback or to requiring that 
the timescale to approach equilibrium be larger than approx- 
imately twice the cloud consumption timescale. We find that 
the amplitude of the oscillations is sensitive to the feedback 
strength and to a lesser extent on the ratio of the delay time 
to the consumption timescale. 

We focus on the molecular or self-gravitating cloud den- 
sity in a galaxy as the most likely variable for the DDE. 
This allows recycling of gas over long periods of time as gas 
is recycled through clouds much faster than it is depleted 
by star formation. The consumption timescale is set by the 
lifetime of molecular clouds. When feedback delay times are 
longer than this timescale we predict episodic star formation 
events and with a period approximately 4 times the delay 
timescale. 

At the present time, there are no compelling constraints 
on either the feedback strength or the delay time, i.e. the 
two key parameters of the DDE model. Thus, it is difficult 
to apply the model rigorously although we suggest avenues 
for further exploration. 

There is more than one candidate for the delay time and 
associated feedback mechanisms, in particular, the timescale 
for supernovae to contribute to turbulence, the timescale for 
spiral density waves to evolve, and the timescale for material 
sent into the halo to return to interact with the disk. We as- 
sociate these three candidate delay mechanisms with possi- 
ble explanations for episodic star formation events in galaxy 
centers (on 10 Myr timescales), the solar neighborhood (on 
Gyr timescales) and dwarf galaxies (on 100 Myr timescales), 
respectively. Using a log normal density distribution we es- 
timate that feedback is likely to be strong enough that the 
second condition for episodic solutions can be satisfied. 

The approach outlined here is potentially powerful 
framework to interpret and motivate future observations and 
simulations. Similar models might be applied to other ac- 
creting systems with feedback such as cooling fiows. With 
better observationally constrained models we may be able to 
use similar simple dynamical models as recipes to drive sim- 
ulations or interpret statistics of astrophysical objects that 
exhibit episodic accretion. 

Lacking currently are simulations and observational 
programs that constrain the timescales and strengths of pos- 
sible feedback mechanisms and their functional form. In view 
of this uncertainty, we adopted an exponential function for 
the feedback process, but is this fully justified? Evidence for 
feedback-infiuenced star formation could be sought by prob- 
ing for correlations between turbulence and deviations from 
empirical star formation laws. Other forms for the feedback 
function could be used, such as that of t he Mackey-Glass 
model which can exhibit chaotic behavior ([Glass fc MackevI 
1988). More sophisticated global theories of star formation 
could be developed to better predict the form of the feed- 
back and go beyond the self-limiting equilibrium state mod- 
els. Higher dimensional models could be explored, similar 
to those used to model predator and prey populations. By 
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going to systems with additional variables it should be pos- 
sible to model these systems without delays. The period is 
not strongly dependent on the amplitude of oscillation for 
the simple model explored here, however, this may not be 
true for more complex models. 
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APPENDIX A: APPLICATION OF 
LINEARIZED OSCILLATION THEORY 

We would like to know when the non-linear DDE given in 
equation [7] exhibits oscillating solutions. For a — 1 this 
differential equation is the same as that of the Lasota- 
Wazewska model. In this appendix we search for a more 
general criterion for oscillation that allows non-unity values 
of the index a. This is desirable because star formation laws 
have non-unity values for this index. 

Non-linear DDEs can have oscillating solutions when an 
associated delay linear equation does. The non-linear DDE 



+ ^^Pifi{x{t - Ti)) 







can be associated with the linearized equation 



(Al) 



(A2) 



(jKulenovic fc Ladasl[l987l : iGvori fc Ladadll99lh . Here > 
0, Ti ^ 0, and the functions fi are well behaved continu- 
ou s functions. Given additional conditions on the fu nctions, 
/-. . iKulenovic fc LadasI ll987l): lGvori fc LadasI proved 
that every solution of the non-linear equation oscillates if 
and only if every solution of the associated linearized equa- 
tion does. One condition is the requirement that 

hm ^ = 1. (A3) 

u^O U 

By manipulating equation [7] and requiring the above 
condition, we find an associat ed linearized equat io n tha t 
is similar to t hat used by IKulenovic fc LadasI |l983); 
iGvori fc Ladad (|l99ll ) to establish when solutions oscillate 
for the Lasota-Wazewska model. This associated linearized 
equation is in the form 



hit) + pix(t) + p2x{t - r) = 0. 



(A4) 



A necessary and sufficient condition for the oscillation of all 
solutions of this linear DDE is 



as proved bv lGvori fc Lada3 (|l99ll ) in section 2.2. Once we 
find the coefficients pi and p2 of the associated linearized 
equation, we can use the above oscillation criterion to es- 
tablish when oscillating solutions exist for the original non- 
linear DDE. 

We wish to find an associated linearized equation for 
the differential equation [7] restated here 

±{t) = Ae"'*-"''''^ 



Bx{tY 

with equilibrium solution, 
change of variables 

x{t) = x,+ Cu(t) 

leads to the delay equation 

Bx'i 



(A6) 

given by equation |3] The 
(A7) 



u{t) + - 



c 



cu(t) y 



+- 



Bx° 



C 



= 0.(A8) 



This can be written in the form of the linearized equation 
[A2lwith 

T-> a — 1 

Bax, 
Bx'i 



Pi 

P2 

h{u) 



c 

Xt 

1-e 



- 1 



(A9) 



where the functions /i , /2 satisfy the condition shown in 
equation IA3I The linearized equation is then in the form of 
equation IA4I pi and p2 into equation IA5|) we find that the 
requirement for oscillating solutions is 



Bx"t laBx°- 



~1) 



> 1. 



(AlO) 



(piT + l) ^ , 



(A5) 



This is dimensionally correct and reduces to equation [TO] 
for the oscillation criterion for the Lasota-Wazewska model 
when a = 1, as expected. 
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